!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine free_energy( energy,iout )
  use rbf_type
  use variable, only: nn, nip, beta, gamma_s, gamma_w1, gamma_w2,ry,rz
  use variable, only: chi, normal, xm, a, da, lhat, ao_s, beta,ao_w
  use variable, only: Sbulk, ULdG, Uinitial, step,q0,ix1,ix1_surf,ddz,pi
use patch_module, only : inside_nano,surface_nano
  implicit none
  real(dp),  dimension(8) :: energy
  integer(i4b) i, j, k, l, ii, m ,iout,imc
  real(dp) vi, ai, tr, tr3, A5, A6, r,tr4
  real(dp) rapini, galatola, tr_gal, UU, SS
  real(dp) aux1, aux2, aux3,k0,k1,k2
  real(dp), dimension(3) :: n, v_gal
  real(dp), dimension(5) :: v1_gal, qo_w
  real(dp), dimension(3,3) :: qq, qow
  real(dp), dimension(5,3) :: grad, v2_gal
   real(dp), dimension(:,:) ,allocatable:: da_local1

   real(dp)    , dimension(3,3)   ::  q,pp,qbar,qdeg

   allocate (da_local1(nn*3,5))

  ! Initializing
  energy = 0.d0

  ! Equilibrium parameters
  SS = Sbulk  ! Only used at the boundary 
  UU = ULdG


  ! differential elements

  vi = (8.d0*ry*rz*beta**3)  /nip
  ai = (ry*8.d0*beta**2 )    /nn 

  A5 = 1.d0/3.d0*SS; A6 = SS*SS
  
  aux1 = 1.d0/sqrt(6.d0)
  aux2 = sqrt(1.5d0) 
  aux3 = sqrt(3.d0)



  ! Internal Energy
  do i = 1, nip
     ii = nn + i
      if(.not.surface_nano(ii).and..not.(inside_nano(ii)))then
     ! Gradients
     grad(1,:) = da(3*ii - 2:3*ii, 1)         ! 1
     grad(2,:) = da(3*ii - 2:3*ii, 2)         ! 2
     grad(3,:) = da(3*ii - 2:3*ii, 3)         ! 3
     grad(4,:) = da(3*ii - 2:3*ii, 4)         ! 4
     grad(5,:) = da(3*ii - 2:3*ii, 5)         ! 5
     


     ! Traces

     

     tr = a(1,ii)*a(1,ii) + a(2,ii)*a(2,ii) + a(3,ii)*a(3,ii) + a(4,ii)*a(4,ii) + a(5,ii)*a(5,ii)



     tr3 = a(1,ii)*a(1,ii)*a(1,ii)*aux1 - aux2*a(1,ii)*a(2,ii)*a(2,ii)  + 0.5d0*aux2*a(1,ii)*a(3,ii)*a(3,ii) &
          + 3.d0*a(3,ii)*a(3,ii)*a(4,ii)/2.d0/sqrt(2.0d0) - aux2*a(1,ii)*a(4,ii)*a(4,ii) + 3.d0*a(2,ii)*a(3,ii)*a(5,ii)/sqrt(2.d0) & 
          + 0.5d0*aux2*a(1,ii)*a(5,ii)*a(5,ii) - 3.d0*a(4,ii)*a(5,ii)*a(5,ii)/2.d0/sqrt(2.d0)


     ! Landau de Gennes
     energy (1)  = energy(1) +  0.5d0*( 1.d0 - UU/3.d0 )*tr - UU/3.d0*tr3 + UU/4.d0*tr*tr
    


!    write(650,*)ii,vi*0.5d0*( 1.d0 - UU/3.d0 )*tr - UU/3.d0*tr3 + UU/4.d0*tr*tr

     
     ! L1 - Elastic
     energy (2)  = energy(2) + 0.5d0*( & 
          grad(1,1)* grad(1,1) + grad(1,2)*grad(1,2)  + grad(1,3)*grad(1,3)  + grad(2,1)*grad(2,1)  + grad(2,2)*grad(2,2)  + &
          grad(2,3)* grad(2,3) + grad(3,1)*grad(3,1)  + grad(3,2)*grad(3,2)  + grad(3,3)*grad(3,3)  + grad(4,1)*grad(4,1)   + & 
          grad(4,2)* grad(4,2) + grad(4,3)*grad(4,3)  + grad(5,1)*grad(5,1)  + grad(5,2)*grad(5,2)  + grad(5,3)*grad(5,3)   )
 

 


     ! L2 - Elastic
     energy (3)  = energy(3)  + 0.5d0*lhat(2)*( & 
          grad(1,1)**2 + grad(1,2)**2 + 4.d0*grad(1,3)**2 - 2.d0*aux3*grad(1,1)*( grad(2,2) + grad(3,3) + grad(4,1) )  + &
          4.d0*aux3*grad(1,3)*( grad(3,1) + grad(5,2) ) + 3.d0*( grad(3,1) + grad(5,2) )**2 - &
          2.d0*aux3*grad(1,2)*( grad(2,1) - grad(4,2) + grad(5,3) ) + &
          3.d0*( ( grad(2,2) + grad(3,3) + grad(4,1) )**2 + ( grad(2,1) - grad(4,2) + grad(5,3) )**2 ) )/6.d0 


     ! L3 - Elastic
!     ener_L3 =  0.5d0*lhat(3)*( &
!          - ( aux3*( grad(1,1)**2 + grad(1,2)**2 - 2.d0*grad(1,3)**2 + grad(2,1)**2 + grad(2,2)**2 - 2.d0*grad(2,3)**2 + &
!          grad(3,1)**2 + grad(3,2)**2 - 2.d0*grad(3,3)**2 + grad(4,1)**2 + grad(4,2)**2 - 2.d0*grad(4,3)**2 + grad(5,1)**2 + & 
!          grad(5,2)**2 - 2.d0*grad(5,3)**2 )*a(1,ii) ) + &
!          3.d0*( 2.d0*( grad(1,1)*grad(1,2) + grad(2,1)*grad(2,2) + grad(3,1)*grad(3,2) + grad(4,1)*grad(4,2) + grad(5,1)*grad(5,2) )*a(2,ii) + &
!          2.d0*( grad(1,1)*grad(1,3) + grad(2,1)*grad(2,3) + grad(3,1)*grad(3,3) + grad(4,1)*grad(4,3) + grad(5,1)*grad(5,3) )*a(3,ii) + &
!          ( grad(1,1)**2 - grad(1,2)**2 + grad(2,1)**2 - grad(2,2)**2 + grad(3,1)**2 - grad(3,2)**2 + grad(4,1)**2 - grad(4,2)**2 + & 
!          grad(5,1)**2 - grad(5,2)**2 )*a(4,ii) + 2.d0*( grad(1,2)*grad(1,3) + grad(2,2)*grad(2,3) + grad(3,2)*grad(3,3) + grad(4,2)*grad(4,3) + &
!          grad(5,2)*grad(5,3) )*a(5,ii) ) )/( 3.d0*sqrt2 )

     energy (4)  = energy(4) + 0.d0

     ! L4 - Elastic
     energy (5)  = energy(5) +  0.5d0*lhat(4)*( &
          grad(1,1)**2 + grad(1,2)**2 + 4.d0*grad(1,3)**2 - 2.d0*aux3*grad(1,1)*( grad(2,2) - 2*grad(3,3) + grad(4,1) ) - &
          2.d0*aux3*grad(1,3)*( grad(3,1) + grad(5,2) ) + 3.d0*( grad(3,1)**2 + 2.d0*grad(3,1)*grad(4,3) - 2.d0*grad(4,3)*grad(5,2) + grad(5,2)**2 ) + &
          2.d0*aux3*grad(1,2)*( -grad(2,1) + grad(4,2) + 2.d0*grad(5,3) ) + 3.d0*( 2.d0*grad(2,3)*grad(3,2) + grad(3,3)**2 + &
          ( grad(2,2) - grad(4,1) )**2 + ( grad(2,1) + grad(4,2) )**2 + 2.d0*( grad(2,3) + grad(3,2) )*grad(5,1) + grad(5,3)**2) )/6.d0     

     energy (7) = energy(7) -  0.5d0*lhat(5)*( (a(5,ii)*(Sqrt(3.d0)*grad(1,1) + grad(2,2) - grad(3,3) + grad(4,1)) + &
          a(4,ii)*(2.d0*grad(2,3) - grad(3,2) - grad(5,1)) + &
          Sqrt(3.d0)*a(1,ii)*(grad(3,2) - grad(5,1)) + &
          a(2,ii)*(grad(3,1) - 2.d0*grad(4,3) - grad(5,2)) + &
          a(3,ii)*(-(Sqrt(3.d0)*grad(1,2)) - grad(2,1) + grad(4,2) + grad(5,3)))/2.d0)

  endif
  end do

360 format (I8,5f20.12)
  ! Surface Energy

 ! Surface Energy
  do i = 1, nn
     if(.not.surface_nano(i).and..not.(inside_nano(i)))then
     ! Normal vector
     n   = normal(3*i - 2:3*i)
     
     ! Galatola vectors
     v1_gal(1) = aux2*( 2.d0*n(3)*n(3)/3.d0 - n(1)*n(1)/3.d0 - n(2)*n(2)/3.d0 )
     v1_gal(2) = sqrt(2.d0)*n(1)*n(2)
     v1_gal(3) = sqrt(2.d0)*n(1)*n(3)
     v1_gal(4) = ( n(1)*n(1) - n(2)*n(2) )/sqrt(2.d0)
     v1_gal(5) = sqrt(2.0d0)*n(2)*n(3)

     v2_gal = 0.d0
     v2_gal(1,1) = -n(1)/3.d0   ; v2_gal(1,2) = -n(2)/3.d0 ; v2_gal(1,3) = 2.d0*n(3)/3.d0 ; v2_gal(1,:) = v2_gal(1,:)*aux2
     v2_gal(2,1) = n(2) ;         v2_gal(2,2) = n(1) ;                                      v2_gal(2,:) = v2_gal(2,:)*sqrt(2.d0)/2.d0
     v2_gal(3,1) = n(3) ;                                    v2_gal(3,3) = n(1) ;           v2_gal(3,:) = v2_gal(3,:)*sqrt(2.d0)/2.d0
     v2_gal(4,1) = n(1) ;         v2_gal(4,2) = -n(2) ;                                     v2_gal(4,:) = v2_gal(4,:)/sqrt(2.d0)
                                  v2_gal(5,2) = n(3) ;       v2_gal(5,3) = n(2) ;           v2_gal(5,:) = v2_gal(5,:)*sqrt(2.d0)/2.d0

     rapini = 0.d0; v_gal = 0.d0; tr_gal = 0.d0

     do j = 1, 5
        rapini = rapini + ( a(j,i) - ao_s(j,i) )**2
        tr_gal = tr_gal + a(j,i)*v1_gal(j)
        v_gal  = v_gal  + a(j,i)*v2_gal(j,:)
     end do

     tr_gal = A5 - tr_gal 
     r        = dot_product( n, v_gal )
     galatola = 2.d0*dot_product( v_gal, v_gal ) + 2.d0*r*r + 4.d0*r*tr_gal + tr_gal*tr_gal
     
     energy(6) = energy(6) + 0.5d0*gamma_w1*( 1.d0 - chi(i) )*galatola + 0.5d0*gamma_s*chi(i)*rapini
     endif
  end do

!  write(*,*)energy(6)

!  energy(6)=0.d0

 ! Surface Energy
  do i = 1, nn
     
     ! Normal vector
     n   = normal(3*i - 2:3*i)
      call inv_transformation_tensor( a(:,i), q (:,:))



     ! tensor
     qq(1,1) = q(1,1) ; qq(1,2) = q(1,2) ; qq(1,3) = q(1,3)
     qq(2,1) = q(2,1) ; qq(2,2) = q(2,2) ; qq(2,3) = q(2,3)
     qq(3,1) = q(3,1) ; qq(3,2) = q(3,2) ; qq(3,3) = q(3,3)
     
     
     ! Qbar = Q + 1/3*S
     qbar(1,1) = q(1,1) + A5; qbar(1,2) = q(1,2);      qbar(1,3) = q(1,3) 
     qbar(2,1) = q(2,1);      qbar(2,2) = q(2,2) + A5; qbar(2,3) = q(2,3)
     qbar(3,1) = q(3,1);      qbar(3,2) = q(3,2);      qbar(3,3) = q(3,3) + A5
     
     ! Proyection
     pp(1,1) = 1.d0 - n(1)*n(1); pp(1,2) =      - n(1)*n(2); pp(1,3) =      - n(1)*n(3)
     pp(2,1) =      - n(2)*n(1); pp(2,2) = 1.d0 - n(2)*n(2); pp(2,3) =      - n(2)*n(3)
     pp(3,1) =      - n(3)*n(1); pp(3,2) =      - n(3)*n(2); pp(3,3) = 1.d0 - n(3)*n(3)

    tr4 = 0.d0
     do j = 1, 3
        do k = 1, 3
           qdeg(j,k) = 0.d0
           do l = 1, 3
              do m = 1, 3
                 qdeg(j,k) = qdeg(j,k)+ pp(j,l)*qbar(l,m)*pp(m,k)
              end do
           end do

           tr4 = tr4 + ( qbar(j,k) - qdeg(j,k) )**2
        end do
     end do
     
     
    
!    energy(6)   = energy(6) +  0.5d0*gamma_w1*tr4  

     
  end do

!      write(*,*)energy(6),'tensor'
!stop
  
  energy(1:5) = energy(1:5)*vi 
  energy (7)  =  energy(7)*vi

  if(nn.ne.0)then 
  energy(6)   = energy(6)*ai
  else
  energy(6) = 0.d0
  endif

  energy(8)   = energy(1) + energy(2) + energy(3) + energy(4) + energy(5) + energy(6) +energy(7)
end subroutine free_energy
